{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f2a808a9",
   "metadata": {},
   "source": [
    "# Causal Inference Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb72ccd0",
   "metadata": {},
   "source": [
    "## Simpson's paradox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "97133256",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pgmpy.models import BayesianNetwork\n",
    "from pgmpy.inference import VariableElimination\n",
    "from pgmpy.factors.discrete import TabularCPD\n",
    "from pgmpy.inference import CausalInference"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae5e7e98",
   "metadata": {},
   "source": [
    "### Model Definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "20fbce7e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._axes.Axes at 0x7f1761cbb700>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKcAAABuCAYAAABC8QOQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAPqklEQVR4nO2de3RVVX7Hv799EkhCQniEMQGigBAcggKCrMJUQRgeo0ApLmtTYo2Ll2EsS5bj4jVZpgYY22kjEtrAIJp2cBDbyqwpD2EAh8XLJeXRFYgYrLkhrARE8s5NbpLz+/WPe68ymIR7Tu4d7mN//srj/h7J/eTsc/bZ2YdEBBpNMKLudQMaTWdoOTVBi5ZTE7RoOTVBi5ZTE7RoOTVBi5ZTE7RoOTVBi5ZTE7RoOTVBS9S9biDUIaJEAH8J4M8AxANoA1ABYI+InL+XvYU6pO+tW4eIUgD8k1JqATPHREVFcWJiIvfo0QPMjIaGBjidzigAopS6wsw/F5H/uNd9hxpaTgsQUS+l1D5mnpKYmNg+b968qNmzZ6N3797fey0z49KlS9i1axdfvHiRiKiJmReKyO/uQeshiZbTR4joWSLaGR8fb6xbt85IT0/3Oba1tRWbN2+WY8eOERH9XkTmiogrgO2GBVpOHyCiVQDenDFjhrz88suklL3ryJKSEuTm5rLL5fqamYeLSJN/Ow0vtJx3gYheAfDWSy+9hKeffrrb+ZxOJ5YuXWo2NDTcZOYh+gjaOVrOLiCiUQAuvvjii7RgwQK/5XU6ncjKymKXy3XENM2ZfkscZuh5zk4gIlJKHRkxYgT7U0wAiIuLQ25urmLmGUQ0z6/JwwgtZ+esJaL73njjDSMQyUeNGoUpU6aIUur9QOQPB7ScnaCUenX69OkUHx8fsBorVqwgEelFRM8GrEgIo+XsACJ6kpn7vvDCCwGt06NHD4wePVqUUusDWihE0XJ2zE8HDhxodjS57m8yMjIUM48IeKEQRN9b7wDDMB4bOXKkpXPN6upqvPPOO6iqqoLT6UT//v2xcePGu8Z5JvOJiMbpe/F/jJazA5g5Zfz48ZZi8vPzMWvWLDz++OMAAIfD4VOcUgpxcXHtTqdzPgAt523oYb0DRCRqyJAhPr/eNE0UFxdj9OjR337NSrznoivV54AIQcvZMRQV5fugYhgGxo4dixUrVmDLli0oKSmxVMwwDADoYa3F8EfL2THS2NhoKSA3NxerV69Gr169kJubi9OnT/sc63K5AMBawQhAn3N2gFKq8fz58wkjR470OYaIkJ6ejvT0dDQ2NsLhcGDSpEk+xdbX1ysAvtscIegjZweIyBfFxcU+Lzo4d+4c2traAAC1tbUoKSnBuHHjfIqtr69He3u7AvBbO72GM/rI2QEicuTKlSvjAPg0nXTy5EkUFhYiJiYG0dHRyMzMxEMPPeRTrY8//hhKqWbTNOu603M4olcldQARJQCoy8nJoYkTJwa0VmZmZntdXd1uEckMaKEQRA/rHSAiDUR0uqioyAxknUuXLqGuri4KwKuBrBOqaDk7QUSWVlRUqE8++SQg+ZkZ69evN5VSfxCRGwEpEuJoOTtBRC4B2Lp582Z2Op1+z19QUCBNTU3tzPyU35OHCVrOrvkpM9/Izs42W1tb/ZZ03759OHz4MInIQhFp9lviMEPL2QUiIsz8bzU1NcaSJUtMfxxBP/roI2zdutX76b5uJwxjtJxdQES5AFaLyJu1tbVfL1y4kO2egzY2NmLlypXme++9JwDWeL7cTEQxfmo37NBTSZ3gEfN1AH8vIrlERAC2AMhOTU3lrKwsw5dppvr6ehQVFeHo0aMiIteZ+cciUuLZxqbW87JYEWkJ0I8Ssmg5O+BOMe/4XjoR/UpEJsXGxvLw4cPVww8/TGPHjkXv3r3R1taG8vJynD17Fl988YVZWVlpKKWqmTlfRDbckUsL2gVazjvoSsw7XtcLwM+JaDoRjWTmBAAEQIioXSlVaZrm/wAoEJFjXeTRgnaClvM2fBUzAHW1oB2gL4g83CsxAUBE6gD08XyqL5I8aDlxb8X0ogX9PhEvZzCI6UUL+sdEtJzBJKYXLeh3RKycwSimFy2om4iUM5jF9KIFjUA5Q0FML5EuaETJGUpieolkQSNGzlAU00ukChoRcoaymF4iUdCwlzMcxPQSaYKGtZzhJKaXSBI0bOUMRzG9RIqgYSlnOIvpJRIEDTs5I0FML+EuaFjJGUlieglnQcNGzkgU00u4ChoWckaymF7CUdCQl1OL+R3hJmhIy6nF/D7hJGjIyqnF7JxwETQk5dRi3p1wEDTk5NRi+k6oCxpScmoxrRPKggb9pgpENBrAQACTocW0zR0bN6QCmAbg1xLEAgT1AwuIKB7AEQA/8HxJi2kTEakjoj5wC1oBoB1ANIAd97CtLgmYnEQ0AEAmgGmGYYwB0MvzrWbTNC8BOArgNyJS0UWaTQASbvv8VCB6jRQ8gn4JYCjc7/3bRHSos/fAs7PeXwCYo5SaREQ/gPsJI+3MfE1ETgL4SEQCsje534d1IpqvlNrIzD/s2bOnmZycjLS0NKNfv34A3M/puXLlCldWVkpLS4uhlPqKmXNF5Nd35JkG4L8BxAFoAVAH4McictGvDUcYRLQdwN/A/Xtth/vhXFNuH96JaBCATUQ0n4iMvn37mg8++GBUamoqoqOj0dLSAofDIQ6Hw6yrqzOIqIWZ/x3AKs85rn969ZecRDRIKXVURNJGjRrFixcvVsOHD+8ypqKiAtu3b+cLFy4QEVUy8zQRKfUM5/8H93DeDOA3AF4REf0IPj9ARDPg/p0mABAAfyci73iOlG8BWNGvXz9zwYIFUXPnzoVSnV83t7a2YteuXThw4IDZ1NQEAD8TkU1+6dMfchLRSwD+JSUlRTZs2GAMGDDAUnxdXR1ycnLMsrIyBSAH7mFnEYAbADICNWxEMkTUG0ABgL+Fe2T6kVJqHxHdt3z5cpo5c6blnDt37sSHH34oRHSJmSeIiKtbPXZXTiJaDeAXzz33HDIzu/ecp71792Lbtm2A+6/5XeijZcDxHEV3EVGf5ORk5OfnG55HbNvi+vXrWLlypel0OquZeaiINNnN1a15TiJaBOAX2dnZ3RYTAObMmYM1a9YA7k1Yy7SYfxI+JaKY+++/H1u3bu2WmACQnJyMHTt2GPHx8f2UUpc9pwq2sC0nEQ0gom3PPPMMnnrKf4/SmTx5MpYtWwYAeUSU5rfEmg5RSh1MSEiI2bRpk9HVuaUV4uLisG3bNkMpNRDANtu92Q5U6mhSUhKysrLspuiUOXPmYOjQoayUOur35JpvIaKFzDwpLy/PiIry76xifHw8Vq5cqQAsJqLxdnLYkpOIHmfm0Rs2bPDpqbp2yMvLM0RkIBE9H6gakY5S6l+nTZsmw4YNC0j+J554AiNGjBCl1Pt24m3JqZT652HDhpkpKSl2wn0iMTERY8aMEaVUbsCKRDBE9BNm7r1s2TLb54S+sHz5csXMIz1zp5awLCcRxTLzhKysrIAdNb0sXbpUMfMwIkoNdK1IQym1MS0tjePi4gJaZ/jw4ejbt68J4B+txto5cs4zDEPGjRtnI9QaqampiImJMQH8dcCLRR4/nDFjxp9kVdrEiRMNwzCmWo2zcxY8u2/fvgwfxT548CD27XM/4tHhcGDIkCEAgEceeQSLFy++a3xKSgrKysqmAfiljV41HUBECQB6Tpo0yVLciRMnsGfPHrS2tsLlcuHJJ59ERkbGXeMmTJiAQ4cO3We1T8tyKqUeHTp0qM9xs2bNwqxZs3Dr1i289tpr2Lx5s6V6aWlpxtWrVx+22mek47lC3g+gFMAfAJwBcBZAJYCZhmFwYmKiz0fOI0eOYO/evVi3bh2SkpLQ3NyMAwcO+BT76KOPQkQMIuovIrd8rWlZTiJKSEhIuPsL76C8vBwPPPCA5bg+ffpARJKJKGjXHQYxAvf6hEkAmgD0BOAC8LVh+H7J4HQ6sWPHDuTn5yMpKQkAEBsbiwULFvgU36NHD++H9wEInJxw/8CWcTgctuQUERARA6i2UzeCGQD3e0UAvA+LNwDUwL2m0+f38fTp00hLS0NycnJ3e2IrL7Ysp4jU19VZXxVVXl4OOxdR1dXVAHBdRO63HBzBeP6DYDeAYgDH4R7S/1dEmoloXnt7+x5fc129ehXdmQttafn2aYk3rMRZlpOZz5SVlY22GlteXo758+dbLYfS0lLTNM0LlgMjHM+61/ROvv17Zla3bt1C//7975orJiYGLpf9BUafffYZiKidmWusxNmZSthfW1traY6TmVFZWYnBgwdbLnbjxg3AvWpe4ydEpFkp1XLqlG//WDB+/HicPHkSNTVut9ra2nDw4EGf6509exZKqSqrfdo559zHzDh16hQmT57sU0BVVRWSkpIQHR1tqVBpaSlcLpcBwNbtL03niMiFQ4cOPTZ37ty7HmjS0tKQkZGB119/HcwM0zQxdepUn2udOXPGNE3zkNUeba3nVEodGzRo0I8KCwsDepdo1apVfPny5cumaXY2PGlsQkR/DuD4zp07kZiYGLA6xcXFWLt2LQAkWZlGAmzeWxeRV65du2Z89dVXdsJ94ubNm/j8888VM68JWJEIRkROKKW+2bJlS0Cn6AoLC02l1AWrYgL25TyvlDqVk5NjMluaHfCZtWvXmkRUKiK/C0gBDZg569NPP6Xi4uKA5N+7dy8qKioUMz9rJ972vVVmnt3Q0MAFBQV+/8vbuXMnrl+/Tsw81d+5Nd8hIvuI6OO8vDy+bbrHL9y8eRPbt28XAP8gIl/ayWFbThFpEJG/Onz4MN5/33/XK/v378fu3bsBIFtELF/haawhIvNcLldNdna22draevcAH6ipqcHy5csZQImI2D4t69aqFBH5LYDFH3zwgbz99tvS3SG+qKgIhYWFALBGRH7VrWQanxCRNmYeWV1d3bBo0SLTM3Vnm9LSUixZsoRbW1vLmHlMd3L561+D5xPR7oSEBCMvL8+wejehqqoK69atM7/55huIyDIRCdotUsIVIkpQSn0mIiMzMjLIl9VGt8PMKCgokMOHD5NS6ggzz+juPkz+3FShl1LqEDNPHjx4sPn8888bd5sHPXfuHIqKisyysjJDKVXMzNNF5KZfGtLYgoh+BuDN2NhYzJw508jMzERMTOcb09XW1uLdd9+V48ePi2mabSKyUET+yy+9BGA7mnFE9JaIPKGUQp8+fcwhQ4YYiYmJRESor6+Hw+For6mpUaZpEhGdEZFXReSEXxvR2IaIYgGsV0otYeaE2NjY9sGDB6uBAweqqKgouFwuXLt2zayqqoLL5TKUUjeZOR/AL0XE9FsfgdoBj4iiAfwEwNNKqcc8O0yQiDQy81kABwDsFZHmgDSg8QtElA7gGQBTDMNIhXtnulbTNL8E8AmA/xQRR0BqB/H2jJoIJ6R2NtZEFlpOTdCi5dQELVpOTdCi5dQELVpOTdCi5dQELVpOTdCi5dQELVpOTdCi5dQELf8Pm7QD1eN+aIgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 153.071x96.378 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "simp_model = BayesianNetwork([(\"S\", \"T\"), (\"T\", \"C\"), (\"S\", \"C\")])\n",
    "simp_model.to_daft(node_pos={\"T\": (0, 0), \"C\": (2, 0), \"S\": (1, 1)}).render()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0f3d089a",
   "metadata": {},
   "outputs": [],
   "source": [
    "cpd_s = TabularCPD(\n",
    "    variable=\"S\", variable_card=2, values=[[0.5], [0.5]], state_names={\"S\": [\"m\", \"f\"]}\n",
    ")\n",
    "cpd_t = TabularCPD(\n",
    "    variable=\"T\",\n",
    "    variable_card=2,\n",
    "    values=[[0.25, 0.75], [0.75, 0.25]],\n",
    "    evidence=[\"S\"],\n",
    "    evidence_card=[2],\n",
    "    state_names={\"S\": [\"m\", \"f\"], \"T\": [0, 1]},\n",
    ")\n",
    "cpd_c = TabularCPD(\n",
    "    variable=\"C\",\n",
    "    variable_card=2,\n",
    "    values=[[0.3, 0.4, 0.7, 0.8], [0.7, 0.6, 0.3, 0.2]],\n",
    "    evidence=[\"S\", \"T\"],\n",
    "    evidence_card=[2, 2],\n",
    "    state_names={\"S\": [\"m\", \"f\"], \"T\": [0, 1], \"C\": [0, 1]},\n",
    ")\n",
    "\n",
    "simp_model.add_cpds(cpd_s, cpd_t, cpd_c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84d92614",
   "metadata": {},
   "source": [
    "### Inference conditioning on T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fc163a22",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6890d595b5cd4adfa75861277cfd474c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6952911a320f4ddc9295692786078edb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------+----------+\n",
      "| C    |   phi(C) |\n",
      "+======+==========+\n",
      "| C(0) |   0.5000 |\n",
      "+------+----------+\n",
      "| C(1) |   0.5000 |\n",
      "+------+----------+\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b186d036897e4ace8fdc23be8ee153d4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fa69a019d7a14821915d56dedcfd956b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------+----------+\n",
      "| C    |   phi(C) |\n",
      "+======+==========+\n",
      "| C(0) |   0.6000 |\n",
      "+------+----------+\n",
      "| C(1) |   0.4000 |\n",
      "+------+----------+\n"
     ]
    }
   ],
   "source": [
    "# Non adjusted inference\n",
    "infer_non_adjust = VariableElimination(simp_model)\n",
    "print(infer_non_adjust.query(variables=[\"C\"], evidence={\"T\": 1}))\n",
    "print(infer_non_adjust.query(variables=[\"C\"], evidence={\"T\": 0}))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af285fcc",
   "metadata": {},
   "source": [
    "### Inference with do-operation on T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9b2d4fa7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "65e8e0dd710b4ca49493d37697e18779",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------+----------+\n",
      "| C    |   phi(C) |\n",
      "+======+==========+\n",
      "| C(0) |   0.6000 |\n",
      "+------+----------+\n",
      "| C(1) |   0.4000 |\n",
      "+------+----------+\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ce1703d2795b4be29fbccd75e848e491",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------+----------+\n",
      "| C    |   phi(C) |\n",
      "+======+==========+\n",
      "| C(0) |   0.5000 |\n",
      "+------+----------+\n",
      "| C(1) |   0.5000 |\n",
      "+------+----------+\n"
     ]
    }
   ],
   "source": [
    "infer_adjusted = CausalInference(simp_model)\n",
    "print(infer_adjusted.query(variables=[\"C\"], do={\"T\": 1}))\n",
    "print(infer_adjusted.query(variables=[\"C\"], do={\"T\": 0}))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64b87d7b",
   "metadata": {},
   "source": [
    "## Specifying adjustment sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8b924b5d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._axes.Axes at 0x7f1760a8ce20>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKcAAACnCAYAAAB0FkzsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUdklEQVR4nO3de3RV9ZUH8O/+nTwvhBAeNkAiTyk+EE0pQgZQM8hAMWBlqdCGJZXXENERfIBIZhwzFB8VUDpJKY/F4iHRNerEIiCTgUFsRFIQpMoiSEgJmCBKCAl5cc/e80cC4rTAPefem55O9mctFjcr9/fbO+bLef7OlUQESnmR+Ws3oNSVaDiVZ2k4lWdpOJVnaTiVZ2k4lWdpOJVnaTiVZ2k4lWdpOJVnRbREESJKANAFAAGoEJFvW6KuCh0iigTQHUAbANUAjkmY732HZctJRB2JaFVERMRxY8wFAGcAfA7gjwC+Mcb4IyIiThLRBiLqFo4eVHCoyTTLsvZblnUeQCOAIwD2AzgKgC3LqjLGfExEY8PSQyjDT0R9jTFvMvNtcXFx9qBBg6yBAwdi4MCBiImJAQA0Njbi008/RVFREfbs2WNXVlZaxpjDzPxzEdkbsmaUK0REAHKJ6BdEFNm/f3+kpKRQamoqEhMTL72vtLQUH3/8Mfbu3cvFxcVEROeZ+UURWRiyXkIVTiJ6EcAzycnJ8uijj5qbb745oHFffvklcnJy+MiRIwRgFYDp4d5dqL+MiAYbY7ZYltVuwoQJ5v7770dExLWP/Gpra7F+/Xq8//77AuA4M98lIqVB9xNsDojIMsZ8BuDG6dOn05gxY1zN8+GHH2LJkiXMzF8xcz8ROR9UY8oRIpoD4NVbb72Vs7KyzMU9nROVlZVYsGCBffz4cQLwUxF5L6iegglnczAPR0VF9czNzTWdOnUKphfU1NRgxowZdk1NzRlmThaRhqAmVAEhoqcBvDx16lSMGzcu6PlycnJky5YtAHBfMAENKpyWZW2PjIy8c+XKlaZ9+/au57lcbW0tpkyZYtfW1n5h2/atIZlUXRERjQTwwfTp05Genh6yeZsDKgB6u93Fuz5bJ6IHmPnuF198MWTBBACfz4clS5ZYInILET0RsonVnyGiSCJ6Z8iQIRLKYAJAZmYmJScnizFmh9s5XIWTiCwiWjdixAjp06eP29pXlJiYiAcffJCI6FdEFB/yAgoAQEQbo6OjY5955hkKx/wLFy60AHQnoufcjHe75ZxrWVbkY489FpYfCgAyMjLg8/kAYFG4arRmRBQtIj/NzMw0gZyRu5GQkIAxY8aQMWaum/GuwmmMmT1s2DBjTHjvfo4ePdoyxjwc1iKt1z9HRUXJ3XffHdYiGRkZEJE4InJ83OA4XUTUg5k7PfLII06HOjZx4kSISCwRhfe/YCtkjHlk2LBhVrjr+Hw+3HDDDUxE852OdbPpeyA6OtoO9CRo6dKlePzxxy/9mTRpEiZOnBjQ2KioKLRr184GcJ+LPtVVMPN1d95557XegwkTJlz6uqSkBOnp6Th58iQAoK6uDg8//DDq6uquOk9KSoohopuc9ujmYOPurl27BvzmJ5544tLriooKzJs3D7NmzQp4fM+ePa0DBw4MddKgajppBXAngC9EpOL/fK8vANO/f/+rzmGMgYhAREBE2Lp1Kzp06IDa2loAwM6dO3HHHXcgNjb2qvOkpqYiLy+vHRGRk7t/jrecxpg+SUlJjncH586dw/PPP4+HHnoId9xxR8DjevToQcaYZKf1FMYB2AKglIjOENEOInqeiMYAGBUZGWkHciIUExODhoYG1NXV4fPPP8eQIUMubSk/+OADjB49+ppz9OzZ8+LL7k5+AMdbTiKKio6OdjSmoaEBL7zwAoYOHRrQD3O56OhoiEhnItL77e5FA7ir+Q8AwBjDgQyMjY1FXV0dPvnkEwwbNgx+vx+1tbU4evQoIiIiLg9eIOKcvNlxOEWk8cKFCwG/37ZtvPzyy0hKSkJGRobTcvD7/SCiMwAC39wqoOk4fSGAegA+AOUAigB8CCCRmZ8OZBKfz4e6ujoUFBRg/vz52LFjB+rq6rB161aMGjXKaU+O1ks4DicznywvL++NAA8JcnNz4ff78dhjjzktBQA4fvw4RKRcRL50NUErRUS5AM4B+ALAfhGpuex7P7pw4cI8Zsa1Lgf6fD4cOHAAHTt2RIcOHRAbG4vKykrs27cPU6dODaiXiopLh7xlTn4GN2frO0+cOBHQLmHjxo04evQo5s2bB8tyd9WipKTEz8yFrga3YiJyXkR+KyIfXR7MZvsAyOHDh685j8/nQ35+/qXDsdjYWGzbtg2DBw9GoId3hYWFMMacF5HAd7lwd7b+dm1t7b/U19fjasuqTp06hTfeeAPXXXcd5s797gZBt27dvvf11TAzKisrLQCbXfSprkBExLKsyl27dnW48cYbr/pen88HZsaAAQMANIXz5MmTeO65wO9I7tu3TwA43vO5WpVkWVZ1enp620A3627l5+dj1apVF0QkWhcghxYRLW/Tps2UvLy8sF6I9/v9GD9+vDDzNBFZ5WSsq/uPzLxy27ZttpuxTrz99tt+EXlXgxkWz5w/f97s378/rEXeeecdiEgjgNVOx7q9Ob6grq6ONm7c6HL4tW3btu3iLn1O2Iq0YiJSZYwpXLx4cdg2MvX19cjLy2MRWe1mA+MqnM2PUMzduHGjnDp1ys0UV1VTU4OcnBwB8JqInAx5AQUAYOYxZ8+elRUrVoRl/hdeeIFt2z4H4FE3410vKxKRXxHR4aeeespubGx0O82fYWbMmTPHFpFT0K1mWIlIlYjMeO+991BUVBTSufPz83Hw4EFi5tFuD8uCWvPGzIPOnTtXPWPGDLu+vj6YqQDg4vVQu6KiopGZU/RYM/xEZDURrczOzpZQBfR3v/sdVq5cCQBPishut/OE4unLjsaYw9HR0QlZWVnXXExwJSUlJViwYIFdU1NTLyI3i8ifgmpMOWKMWSEiU8eOHYtp06a5msPv9+Oll17i3bt3GwBPi8ivgukpJM+tNz+Lki8iowcPHiyzZs2i+PjAnq6ora3Fb37zG9mxYwcZYwqZeZSIVAfdlHKMiB4houXt27enOXPmWLfddlvAY3fs2IGcnBxuaGioE5H7RKQg6H5C/IkfY4wxa5i5U48ePeyRI0daqamp6Nix4/feV1VVhcLCQhQUFHBxcbExxpxj5kwR2RCyZpQrRBRvjNnEzH/Xpk0bTktLs4YOHYp+/fp971an3+/HZ599hp07d+Kjjz6yGxsbiYjeFZGJTu8EXbGXcBzWEdFQInqZiG5j5lhjjERERAjQtBDEtm1jjGkAcIiZnxWRrSFvQgWl+cHCl4wxDzBzAgCKjIy0jTFk2zb8fr8BwMaYr5l5FYB/DVUoL/UQ7nMOIooFcA+A9WhaMjUewAf6iR5/W4jodgD/gKYHDn8BoFBEisNas6VOiInoawCdRSRsT2yq8CKiPgCOtNTvUD88VnmWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeZaGU3mWhlN5loZTeVaLhJOIcgF0bn79VkvUVKFFRMMA7Gp+fYiIeoW7ZkttOVMue/2jFqqpQqsbgHbNr7sDaB/ugi0Vzg8BcPPrT1qopgqtvQCk+XUEgD+Gu2BEuAs0+wRATXO9D1uopgqto/huY1YqIo3hLthSW869aAqmv/m1+hsjIgzgUPOXhS1Rs6XCWYqmXUIbAAdbqKYKvV3Nf/++JYq1SDhFRAB8AeBPIlLfEjVVWOxu/rtF9n7UlJswTEzUGUAGgDTLsgYA6ISmfwzf2Lb9OYDtAN4QkbKwNKCCRkQEYByAe40xQ4ioi4jEE9EZZi4Tkd8DeEdEdoSlfqjDSUT3GWN+ycw3RkdH24mJiejbt6/VoUMHAMDZs2dx5MgR/uqrr6S+vt4yxpQw8/Misi6kjSjXiKgbgKVEdB8RWQkJCXbv3r0jkpOTERkZifr6epSWlkppaaldVVVlEVE9M68FMFdEqkLWR6jCSUTdjDHbRaTvTTfdxFOnTjV9+vS56piysjKsWLGC9+/fT0T0FTOniUhxSBpSjjVvKZcAeLxDhw72/fffH5Geng5jrnz019jYiI0bN2LLli32+fPnAeApEVkakn5CEU4i+kcA/96lSxdZuHCh1blzZ0fjq6qqkJWVZR87dswAyBKRhUE3pRwhoi7GmH1E9IPMzEwaOXKk4znWr1+Pt956S4joc2YeKCINQfUUbDiJaB6ARQ899BAyMjKCmmvTpk1Yvnw5ALwuIv8U1GQqYETUnYgOJSYmRi1evNhq27at67kqKiowe/Zsu7a29gwz9xSR8677CiacRDQFwMqZM2fiJz/5iet5LldYWIhFixYBwALdgoYfEcURUfn1118f8/rrr1tX24UHqra2FtOmTbNramrKmfl6cRky150QUWciWj5+/PiQBRMAUlNTMWPGDADIJqK+IZtY/UXGmA/i4uJili5dGpJgAoDP58Py5cstY0xXAMtd9+Z6oDHbO3XqhMmTJ7ud4oruvfde9OzZk40x20M+ubqEiH7OzEOys7OtiIjQ3slu27YtZs+ebQBMJSJXi31chZOIhjHzLQsXLrTcjA9Edna2JSJdiWhSuGq0dsaYnLS0NOnVKzyr34YPH44bbrhBjDEb3Ix3FU5jzKu9evWyu3Tp4mZ4QOLj4zFgwAAxxjwftiKtGBGNZuZ2M2bMoHDWyczMNMz8w+Zrp444DicRxTLzwMmTJ4dtq3nR9OnTDTP3IqLkcNdqbYwxv+zbty/7fL6w1unTpw8SEhJsAC87HetmyznWsiy5/fbbA3rz5s2bkZOTc+nrdevW4dVXXw1obHJyMmJiYmwAE1z0qa7uxnvuuSeg339ZWRmmTJkC5qYlucyMrKwsbN8e2CnBoEGDLMuy7nLaoJtwjkpISOBrv61JWloaioqKUFNTgz179uAPf/gDZs2aFXCx5kOHNBd9qisgojhmjh4yZEhA709OTkZSUhKKiooAAGvXrkVSUhLS0gL7tQwcOBDM/AOnfTo+RTPGpPTs2TPgcTExMRg+fDjWrVuHvXv3Ijs7G9HR0QHX69u3r3X8+PH+Tvts7ZrPkDcDKAbwPwCK0LSa6CsAIy3L4vj4+IA3TuPGjUN+fj78fj8OHTqEhQsDvwSdkpICEbGIqKOIfBvoOMfhJKK4uLg4R2NGjBiBzMxMLFiwAE5Potq3bw8RSSSi8Cyf+v9NAFwHYAiA8wCiATQA+NqynJ0ypKSkYPXq1Vi7di0WLVoEJ5eeoqKiLr78AYDwhRPfPUcSsLy8PMTHx8O2befFREBEDOCM48GtW2c0/a4IwMVHKiwAlQDOwsXvsV+/fujduzcurjBzIeDDQcBFOEXkXFVV4Kui3n33XVy4cAFz587Fhg0bkJqa6qjemTNnAKBCRK531mnrRkS3AHgTTU8e7ELTLv2AiNQR0Vi/3/+u0znLysowYsQIx73U119aX37KyTjH4WTmomPHjt0SyNgDBw6goKAAr7zyCnw+H+rq6lBSUgInF32Li4tt27b3O+2ztRORPwK4+Qrf/i9mNt9++y06duwY8JzHjx9H9+7dHfeyZ88eEJGfmSudjHNztr757Nmz1zxg+frrr7Fs2TLMmzcPF6+lpaenIz8/31GxU6dOAU2r5lWIiEidMaa+sDDw59ROnz6NNm3aIDY21nG9vXv3whhT7nSc41VJRBQJoOHZZ58lp7top4qLi/Hkk08CwHUicjqsxVoZY8zH3bt3//GyZcvCfjPlZz/7mV1dXb1GRKY6Ged4yykiF4ho17p165yf3Ti0atUqNsZ8ocEMPRF5urS01HJy/uDGwYMHUV1dbQGY63Ssq3vrIvLEiRMnrJKSEjfDA3L69GkcOnTIMPOzYSvSionIR8aYb37961+H9RJdbm6ubYzZ7+T65kVuw/mpMaYwKyvLvnhLK9Tmz59vE1GxiLwXlgIKzDx59+7ddPBgeD5KYNOmTSgrKzPM/ICb8a7XczLzqOrqal62bFnI/+WtX78eFRUVxMx3hXpu9R0ReZ+ItmZnZ/Nll3tC4vTp01ixYoUAeElEvnQzh+twiki1iDxYUFCADRtcLdf7izZv3ow333wTAGaKiOMzPOWMiIxtaGionDlzpt3YGJqPP6qsrERmZiYD+EJEXB+WBbUuX0T+E8DUvLw8ee211yTYXfyaNWuQm5sLAM+KyG+DmkwFREQuMPMPz5w5Uz1lyhS7+dKda8XFxZg2bRo3NjYeY+YBwcwVqkeD7yOiN+Pi4qzs7GzL6crq8vJyPPfcc/Y333wDEZkhIquCbko5QkRxxpg9IvLDiRMn0sSJEx2NZ2YsW7ZMCgoKyBjz38x8j9sH2y71FMIPVWhjjNnGzKlJSUn2pEmTrGtdB923bx/WrFljHzt2zDLGHGTmv9fLRn9dRPQUgBdjY2MxcuRIKyMjAzExMVd8/9mzZ7F69WrZtWuX2LZ9QUR+LiJvh6SXMHwcze1EtEREhhtj0L59e7tHjx5WfHw8ERHOnTuH0tJSf2VlpbFtm4ioSESeFJGPQtqIco2IYgH8mzFmGjPHxcbG+pOSkkzXrl1NREQEGhoacOLECbu8vBwNDQ2WMeY0My8G8IqIhOz6dzg/yCsSwGgAY4wxPyaids31aph5L4AtADaJSF1YGlAhQUQ3AxgP4E7LspIBRAJotG37SwA7APyHiJSGpXa4wqlUsPR/9aI8S8OpPEvDqTxLw6k8S8OpPEvDqTxLw6k8S8OpPEvDqTxLw6k8S8OpPOt/AW2V9Noo7r6bAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 153.071x153.071 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model = BayesianNetwork([(\"X\", \"Y\"), (\"Z\", \"X\"), (\"Z\", \"W\"), (\"W\", \"Y\")])\n",
    "cpd_z = TabularCPD(variable=\"Z\", variable_card=2, values=[[0.2], [0.8]])\n",
    "\n",
    "cpd_x = TabularCPD(\n",
    "    variable=\"X\",\n",
    "    variable_card=2,\n",
    "    values=[[0.1, 0.3], [0.9, 0.7]],\n",
    "    evidence=[\"Z\"],\n",
    "    evidence_card=[2],\n",
    ")\n",
    "\n",
    "cpd_w = TabularCPD(\n",
    "    variable=\"W\",\n",
    "    variable_card=2,\n",
    "    values=[[0.2, 0.9], [0.8, 0.1]],\n",
    "    evidence=[\"Z\"],\n",
    "    evidence_card=[2],\n",
    ")\n",
    "\n",
    "cpd_y = TabularCPD(\n",
    "    variable=\"Y\",\n",
    "    variable_card=2,\n",
    "    values=[[0.3, 0.4, 0.7, 0.8], [0.7, 0.6, 0.3, 0.2]],\n",
    "    evidence=[\"X\", \"W\"],\n",
    "    evidence_card=[2, 2],\n",
    ")\n",
    "\n",
    "model.add_cpds(cpd_z, cpd_x, cpd_w, cpd_y)\n",
    "\n",
    "model.to_daft(node_pos={\"X\": (0, 0), \"Y\": (2, 0), \"Z\": (0, 2), \"W\": (2, 2)}).render()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c0b2e092",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "64df707c71d34aa192f465fffcdd2d1d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------+----------+\n",
      "| Y    |   phi(Y) |\n",
      "+======+==========+\n",
      "| Y(0) |   0.7240 |\n",
      "+------+----------+\n",
      "| Y(1) |   0.2760 |\n",
      "+------+----------+\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "351bed21a88b49d8ac1d06c72a0ce751",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------+----------+\n",
      "| Y    |   phi(Y) |\n",
      "+======+==========+\n",
      "| Y(0) |   0.7240 |\n",
      "+------+----------+\n",
      "| Y(1) |   0.2760 |\n",
      "+------+----------+\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "193a7875ffcb4a80a193e7d1d1d57f76",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/4 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------+----------+\n",
      "| Y    |   phi(Y) |\n",
      "+======+==========+\n",
      "| Y(0) |   0.7240 |\n",
      "+------+----------+\n",
      "| Y(1) |   0.2760 |\n",
      "+------+----------+\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "aabdcce521ca42d8b7180019aac16b66",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/2 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------+----------+\n",
      "| C    |   phi(C) |\n",
      "+======+==========+\n",
      "| C(0) |   0.6000 |\n",
      "+------+----------+\n",
      "| C(1) |   0.4000 |\n",
      "+------+----------+\n"
     ]
    }
   ],
   "source": [
    "# Do operation with a specified adjustment set.\n",
    "infer = CausalInference(model)\n",
    "do_X_W = infer.query([\"Y\"], do={\"X\": 1}, adjustment_set=[\"W\"])\n",
    "print(do_X_W)\n",
    "\n",
    "do_X_Z = infer.query([\"Y\"], do={\"X\": 1}, adjustment_set=[\"Z\"])\n",
    "print(do_X_Z)\n",
    "\n",
    "do_X_WZ = infer.query([\"Y\"], do={\"X\": 1}, adjustment_set=[\"W\", \"Z\"])\n",
    "print(do_X_WZ)\n",
    "\n",
    "infer_simp = CausalInference(simp_model)\n",
    "do_simpson = infer_simp.query([\"C\"], do={\"T\": 1}, adjustment_set=[\"S\"])\n",
    "print(do_simpson)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "98ea8930",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------+----------+\n",
      "| Y    |   phi(Y) |\n",
      "+======+==========+\n",
      "| Y(0) |   0.6200 |\n",
      "+------+----------+\n",
      "| Y(1) |   0.3800 |\n",
      "+------+----------+\n",
      "+------+----------+\n",
      "| Y    |   phi(Y) |\n",
      "+======+==========+\n",
      "| Y(0) |   0.6200 |\n",
      "+------+----------+\n",
      "| Y(1) |   0.3800 |\n",
      "+------+----------+\n",
      "+------+----------+\n",
      "| Y    |   phi(Y) |\n",
      "+======+==========+\n",
      "| Y(0) |   0.6200 |\n",
      "+------+----------+\n",
      "| Y(1) |   0.3800 |\n",
      "+------+----------+\n",
      "+------+----------+\n",
      "| C    |   phi(C) |\n",
      "+======+==========+\n",
      "| C(0) |   0.5500 |\n",
      "+------+----------+\n",
      "| C(1) |   0.4500 |\n",
      "+------+----------+\n"
     ]
    }
   ],
   "source": [
    "# Adjustment without do operation.\n",
    "infer = CausalInference(model)\n",
    "adj_W = infer.query([\"Y\"], adjustment_set=[\"W\"])\n",
    "print(adj_W)\n",
    "\n",
    "adj_Z = infer.query([\"Y\"], adjustment_set=[\"Z\"])\n",
    "print(adj_Z)\n",
    "\n",
    "adj_WZ = infer.query([\"Y\"], adjustment_set=[\"W\", \"Z\"])\n",
    "print(adj_WZ)\n",
    "\n",
    "infer_simp = CausalInference(simp_model)\n",
    "adj_simpson = infer_simp.query([\"C\"], adjustment_set=[\"S\"])\n",
    "print(adj_simpson)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}